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SECTION  1 


INTRODUCTION 

The  inviscid,  two-invariant  associative  cap  model  was  originally  proposed  by  DiMaggio 
and  Sandler  [1,2]  and  an  algorithm  to  implement  the  model  in  stress  analysis  programs  was 
proposed  by  Sandler  and  Rubin  [3],  To  assess  the  predictive  capabilities  of  the  inviscid  cap 
model,  the  extensive  and  well-documented  data  obtained  in  the  experimental  program  at  the 
University  of  Colorado  [4]  has  been  selected.  A  characteristic  of  this  experimental  work  is 
the  exercise  of  truly  three  dimensional  non-conventiona!  stress  paths.  Due  to  the  non- 
conventional  nature  of  the  experimental  data,  standard  fitting  procedures  based  on  the  use  of 
conventional  tests  to  independently  fit  the  cap  surface,  failure  envelope  and  hardening  law 
(see  e.g.  [5.6])  cannot  be  used.  Hence,  to  obtain  values  for  the  cap  parameters,  an  alternative 
constrained  optimization  procedure  which  employs  a  modified  Marquardt-Levenberg  algo¬ 
rithm  and  Armijo  step-size  rule  is  developed.  This  approach  makes  the  fitting  process  com¬ 
pletely  systematic  and  renders  the  optimal  values  of  the  parameters  in  a  least  square  sense. 

In  the  simulations  reported  herein,  six  (6)  tests  are  used  to  fit  the  seven  parameters  of 
the  cap  model,  and  the  resulting  model  is  exercised  to  predict  the  remaining  sixty-one  (61) 
tests.  The  resulting  numerical  predictions  agree  remarkably  well,  both  qualitatively  and  quanti¬ 
tatively,  with  the  experimental  results. 


I 


.*»'*<*  I  IaM  *<l  4.1' 


SECTION  2 

BASIC  FORMULATION  OF  THE  INVISCID  CAP  MODEL 

The  two-invariant,  rate-independent  elastoplastic  associative  cap  model  is  characterized 
by  the  following  constitutive  equations: 

t  =■  te  +  tp 

a  =  a(tc)  (elastic  response) 

ip  =  \  (associative  flow  rule)  ( 1 ) 

09 

<j>(a,n)  <,  0  (yield  condition ) 

where  «,  te,  and  tp  denote  the  total,  elastic  and  plastic  strain  tensors;  a  denotes  the  stress  ten¬ 
sor  and  <t>(a,  k)  =  0  is  the  yield  surface  in  stress  space.  In  addition,  x  is  the  hardening  parame¬ 
ter  which  for  the  cap  model  is  related  to  the  plastic  volume  change  by  a  hardening  law  as 
described  below.  Loading/unloading  conditions  may  be  expressed  in  a  compact  manner  by 
requiring  that 

</>(<r.O<  0.  A>0.  A0(».k)  =  O  (2) 

This  is  the  so-called  Kuhn-Tucker  form  of  unilateral  constraint  conditions.  Note  that  if  <i>  <  0 
then  X  =  0  and  the  process  is  elastic.  On  the  other  hand,  for  loading,  X  >  0  and  <t>  =  0.  In 
this  latter  case.  X  is  determined  by  requiring  that  0  =  0;  the  so-called  consistency  condition 
leads  to  the  classical  elastoplastic  tangent  modulus. 

The  basic  characteristic  of  the  cap  model  is  the  form  of  the  yield  function  4>(o,k)  which 
is  specified  in  terms  of  two  functions  I',  and  l\  .  The  function  /■',  denotes  the  so-called  failure 
envelope  surface  whereas  the  function  /•'  is  referred  to  as  the  hardening  cap.  Functional 
forms  for  Ft.  and  I\  arc  (sec  Fig.  I) 


0(*.»O  s 


\  J id  —  h,.(J i )  S  0 

\  J id  —  F*t  (J  ] . . )  ^  0 


( failure  envelope) 
(cap  surface) 


where  J  \  =  ira  .  J m  =  '/.’s  s  (s  :  stress  deviator)  and 


h\  (J \)  =  t<  -  y  cxp(-d  J |)  +  07 1 


i\(J i.O  =  -=-\  (.V(0-L(k)]:-[7,-Z.(»)]- 

A 


L(k)  =  <k> 


k  il  k  >  0 

0  It  k  <,  0 


ft 


(  <MeAuley  hracket>  ) 
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Finally,  the  hardening  parameter  k  is  related  to  the  plastic  volume  change  =  tr  tp  by 
the  hardening  law 

«?(*>■  W  {  1  -  exp  [-£>*(*)]  >  (5) 

where  X(k)  is  defined  by 

X(k)  =  k  +  R  F6{k)  (6) 

In  the  above  expressions,  a  ,0 ,y  ,Q ,  IV  ,D  ,  and  R  are  material  parameters  which  characterize 
the  two-invariant  cap  model  considered  here. 


SECTION  3 


PARAMETER  ESTIMATION  AND  NUMERICAL  SIMULATIONS 

In  order  to  assess  the  capability  of  the  two-invariant  cap  model  in  predicting  response 
behavior  for  actual  materials  such  as  concrete  and  geomaterials,  model  parameters  need  to  be 
estimated  from  available  experimental  data.  In  this  section,  a  parameter  estimation  pro¬ 
cedure  and  an  assessment  of  the  predictive  capability  of  the  cap  model  are  presented.  This  is 
followed  by  extensive  numerical  simulations  for  the  Colorado  concrete  data. 

3.1.  PARAMETER  ESTIMATION.  MARQUARDT-LEVENBERG  ALGORITHM. 

It  is  characteristic  of  currently  employed  parameter  estimation  procedures  for  the  cap 
model  (see  e.g.  [5,6])  to  fit  separately  the  failure  envelope,  cap  surface,  and  hardening  law 
parameters.  Typically,  asymptotic  failure  points  from  TE,  TC,  SS,  CTC,  CTE,  RTE,  RTC 
and  PLf  arc  used  with  a  least-square  fit  procedure  to  estimate  the  failure  parameters;  whereas 
iso-plastic  volumetric  strain  contours  are  employed  to  estimate  the  cap  shape  parameter  R. 
The  hardening  law  parameters  D  and  W  are  determined  from  HC  tests.J  Although  this  pro¬ 
cedure  provides  a  parameter  fitting  directly  associated  with  the  physical  construction  of  the 
cap  model,  it  has  the  following  two  major  drawbacks:  (a)  a  large  amount  (more  than  20  tests) 
of  conventional  experimental  data  are  required  (e.g.  CTC,  CTE  etc.),  and  (b)  it  is  not  possible 
to  utilize  some  existing  nonconventional  experimental  work:  e.g.,  the  results  from  the 
"Colorado"  experimental  program  [4],  Hence,  a  more  flexible  and  systematic  parameter  esti¬ 
mation  procedure  is  needed.  This  is  the  objective  of  the  following  section. 

Optimization  algorithm  The  basic  idea  of  the  procedure  advocated  here  is  to  regard  the 
optimal  fitting  process  for  given  experimental  data  as  a  least-square  constrained  optimization 
problem.  In  this  context,  the  objective  function  11:  R'-^IR  is  simply  the  sum-of-squares 
error  function  defined  as 

,\ 

n(*)s  2  Ik/ (*.«/)  -  */'ll:  (7) 

i = i 

where 

,V  :  number  of  observation  s 

♦  I  L  stands  lor  iria.xial  extension.  IV  iria.xial  compression.  SS  simple  shear.  C  TC  conventional  iria.xial 
compression.  CTE  conventional  iria.xial  extension.  KTE  reduced  tnaxial  extension.  RTC  reduced  triaxia! 
compression,  and  PL  proportional  loading. 

}  HC  represents  hsdroslaltc  compression  test. 


9 1 :  stress  response  from  constitutive  model  considered 


*]\  observed  stress  response 

parameter  vector  (in  R7  for  cap  model) 

/ :  l'1' data  point 

In  the  following,  this  procedure  will  be  illustrated  using  the  cap  model.  It  is,  however,  gen¬ 
erally  applicable  to  any  constitutive  model.  The  constraints  imposed  on  the  optimization 
problem  emanate  from  physical  restrictions  placed  on  the  cap  parameters.  For  example,  for  a 
physically  meaningful  model  one  should  have 

a  >  0,  y  >  0,  a  >  y,  0  >  0,  0  >  0,  R  >  0,  D  >  0,  W  •>  0.  These  constraints  define  a  feasible 
domain  ZC1R7,  which  is  a  convex  polygon.  The  resulting  constrained  optimization  problem  is 
then  expressed  as 

Find:  min  Ilf*)  subject  to  ♦eS  (8) 

There  exists  a  wide  variety  of  algorithms  for  solving  the  standard  convex  optimization 
problem  (8)  (e.g.,  see  [9]  for  a  review).  The  algorithm  employed  here  is  the  well-known 
Marquardt-Levenberg  algorithm  together  with  the  Armijo  step-size  rule  [7-10],  This  algorithm 
is  essentially  a  hybrid  of  Newton  and  steepest  descent  (gradient)  methods.  It  combines  the 
ability  of  the  steepest  descent  method  to  converge  from  an  initial  guess,  which  may  be  outside 
the  region  of  convergence  of  other  methods,  with  the  asymptotic  quadratic  convergence 
characteristics  of  Newton’s  method  near  the  solution.  The  Marquardt-Levenberg  algorithm 
can  be  summarized  in  the  following  form: 

*/♦■-*,•+  \  h,  (9) 

h,  =  -[H,  +  n,  D,  ]- 1  v,n  (10) 

H,  =  2  Q,r  Q,  (approx.  Hessian)  (11) 

Q,  =  (sensitivity  matrix)  (12) 

o'*  I 

i?,  =  Marquardl  parameter 

D,  =  diagonal  matrix  or  simply  I 


X,  =  argmin  w4  I  II(*,.l)<ri(+l) 

..(llll  I  •  [S 


i  -i  iteration 

For  problems  where  Q ,  may  not  be  easily  constructed  analytically  the  derivatives  are  typically 
computed  by  means  of  forward  differences.  However,  central  differences  provide  greater 
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accuracy  in  the  vicinity  of  the  solution  (minimum);  thus,  central  rather  than  forward 
differences  are  employed  in  computing  Q,  when  the  solution  is  closely  approached. 

In  addition,  to  minimize  the  number  of  function  evaluations  (stress  responses),  a  rank 
one  update  to  the  sensitivity  matrix  is  used  periodically  (similar  to  the  Quasi-Newton  method) 


Q,.i  =  Q,  +  — r1— rrd -Q,  r.,  (14) 

I  r 

where  =  ♦,„!  In  Eq.  (10),  for  a  given  value  of  17, ,  Cholesky  factorization  of 

H,  +  n,  D,  is  employed  to  check  for  positive  definiteness.  If  the  factorization  breaks  down,  i.e. 
H,  +n,  D,  is  not  positive  definite,  then  ?j,  is  increased.  The  algorithm  summarized  above  can 
be  systematically  applied  to  any  set  of  experimental  data  to  obtain  the  optimal  fit  for  the  con¬ 
stitutive  model  under  consideration  in  a  least  square  sense. 

Error  measurement  During  the  optimization  process,  a  root-mean-square  (RMS)  type  of 
error  measurement  is  adopted.  The  optimization  process  is  considered  to  reach  its  optimum 
when  the  RMS  measure  is  minimized.  The  relevant  measures  are  defined  as  follows: 

r  -i  i 


A.v 


U  2 
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( RMS  of  error) 


(15) 
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(RMS  of  observed  responses) 


(16) 


«.v 


(normalized  relative  RMS  error) 


(17) 


Remark  3.1.  It  is  interesting  to  examine  the  sensitivity  of  the  response  under  perturba¬ 
tions  in  cap  model  parameters.  A  finite  difference  sensitivity  matrix  Q  is  defined  in  dimen¬ 
sionless  form: 


Qu 


Act,  /  gy 

A*,  /  *j 


(18) 


where  is  a  stress  component  (i  =  1,...,6)  and  is  a  parameter  component  (j  =  1 . 7), 

respectively.  A  standard  sensitivity  analysis  reveals  that  the  response  of  the  cap  model  is  rela¬ 
tively  insensitive  to  changes  in  the  model  parameters.  By  ordering  the  mode!  parameters 
according  to  relative  sensitivity  in  the  response,  one  obtains  in  decreasing  order  of  sensitivity: 


U  -»  D  —  R  —  ,r— 0—  7— # 


(19) 


In  summary,  one  obtains  the  following  relative  degree  of  sensitivity  (from  large  to  small): 


harden m v  parameters 


eap  parameters 


failure  parameters  □ 


3.2.  PREDICTIVE  CAPABILITIES.  “COLORADO”  CONCRETE  DATA. 

In  this  section,  we  first  examine  the  consistency  of  the  'Colorado  concrete'  data  [4],  next 
we  estimate  the  model  parameters  by  exercising  the  procedure  described  above,  finally  we 
assess  the  predictive  capability  of  the  inviscid  cap  model. 

Colorado  concrete  data.  This  experimental  program  on  concrete  was  performed  at  the 
University  of  Colorado  (1983)  and  is  well-documented  [4],  The  program  consists  of  six  major 
series  of  nonconventional  multiaxial  stress-strain  curves.  The  total  number  of  experiments  is 
67.  The  data  are  characterized  by  the  following  properties:  (a)  characteristic  uniaxial  compres¬ 
sive  strength  f'c~4  ksi,  (b)  mean  pressure  £  8  ksi  (c)  truly  triaxial  states  of  stress  for  con¬ 
crete,  (d)  nonconventional  complicated  stress  paths,  and  (e)  quasi-static  loading. 

The  six  major  series  of  tests  consist  of  the  following: 

(1)  A  series  of  12  cyclic  triaxial  tests,  consisting  of  cyclic  hydrostatic  preloading  to  various 
stress  levels,  followed  by  proportional  deviatoric  stress  cycles  without  reversal  along 
triaxial  compression,  simple  shear,  and  triaxial  extension  paths. 

(2)  A  series  of  8  cyclic  triaxial  tests,  consisting  of  cyclic  hydrostatic  preloading  to  various 
stress  levels,  followed  by  proportional  deviatoric  stress  cycles  with  reversal  along  the 
same  deviatoric  paths  as  in  Series  I . 

(3)  A  series  of  17  tests  consisting  of  hydrostatic  loading,  followed  by  proportional  stress 
deviation,  followed  by  a  circular  stress  path  within  the  deviatoric  plane. 

(4)  A  series  of  22  axisymmetric  triaxial  tests  to  explore  load  path  effects.  In  addition  to  pro¬ 
portional  and  hydrostatic-deviatoric  paths,  this  series  contained  staircase-type  loadings 
to  explore  convergence  to  the  proportional  path,  tests  with  hydrostatic  stress  increments 
with  and  without  hydrostatic  preloading,  and  tests  under  non-proportional  loadings. 

(5)  A  series  of  6  tests  within  the  deviatoric  plane,  as  well  as  a  number  of  other  tests 
specifically  designed  to  check  the  meaning  of  loading  and  unloading. 

(6)  A  series  of  2  tests  of  piecewise-uniaxia!  loadings. 

Assessment  of  data  consistency.  Basically,  the  measures  employed  here  are  the  same  as 
those  discussed  in  the  previous  section.  For  convenience,  these  measures  are  summarized  as 
follows: 
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Here  «/  refers  to  a  strain  measurement  of  test  ’A'.  An  assessment  of  consistency  for  the 
'Colorado'  concrete  data  may  be  obtained  from  the  replicates  of  experiments  available  in  the 
reported  results  [4],  The  present  analysis  generally  indicates  reasonable  consistency  of  the 
data.  However,  some  serious  discrepancies  between  replicates  are  also  observed.  See  Table  1 
below. 


Table  1.  Consistency  of  the  Colorado  concrete  data  [4], 


Tests 

5  % 

Major  Path 

1-1  &  1-10 

13.5 

CTC 

1-4  &  1-7 

31.1 

TC 

1-6  &  1-9 

51.3 

TE 

2-3  &  2-4 

9.6 

SS 

2-7  &  2-8 

13.5 

SS 

3-1  &  3-2 

244.3 

Circular 

3-3  &  3-4 

47.2 

Circular 

3-10  &  3-11 

92.9 

Circular 

4-1  &  4-2 

10.9 

Axisymmetric 

4-6  &  4-7 

54.2 

Axisymmetric 

Model  parameter  estimation  procedure.  The  actual  data  employed  in  the  optimization 
process  based  on  the  Marquardt-Levenberg  algorithm  are  obtained  by  arbitrarily  selecting  one 
test  out  of  each  of  the  six  major  series.  Thus,  a  total  number  of  6  tests  is  used  in  the  actual  fit 
of  the  model.  The  quality  of  the  fitting  is  satisfactory.  Typical  values  of  the  RMS  error 
found  from  back-prediction  using  the  optimal  material  properties  are:  5  =  16%  for  test  1-1 
(CTC),  5  =  8.5%  for  test  2-3  (SS),  5  =  26%  for  test  3-11  (circular),  5=11.5%  for  test  4-11 
(axisym.),  etc..  From  this  optimization  procedure,  we  obtain  the  following  set  of  parameters 
which  best  fits  the  observed  experiments:  a  =  3.86 ksi ,  0  =  .11,  7=1.16fcst,  /?=.44 ksi~\ 

*=4.43,  D  =.0032 ksC\  W  =  A2,  X°  =  \(>ksi. 

Predictive  capability.  After  the  optimal  model  parameters  are  obtained,  the  resulting  cap 
model  is  used  to  predict  the  response  of  every  other  Colorado  test  which  is  not  included  in  the 
optimization  process  (total  number  =  61).  It  is  emphasized  that  the  "prediction"  here  has 
nothing  to  do  with  optimal  fitting,  but  is  obtained  by  exercising  the  cap  model  using  previ¬ 
ously  estimated  parameters.  In  general,  considering  the  experimental  data  scatter,  the 
predicted  response  is  in  good  agreement  with  the  experimental  results.  It  is  noted  that  the 
overall  qualitative  behavior  for  the  Colorado  concrete  data  is  captured.  Values  of  the  RMS 
error  corresponding  to  a  selected  sample  of  simulations  are  summarized  in  Table  2  below. 
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The  overall  RMS  and  standard  deviation  of  error  for  61  tests  are  26.6%  and  14%,  respec¬ 
tively.  A  comparison  between  experimental  and  predicted  stress-strain  curves  is  contained  in 
Figures  2-11. 


Table  2.  Results  of  prediction.  Inviscid  case 


Tests 

5  % 

Major  Path 

1-2 

12.4 

SS 

1-3 

14.1 

TE 

2-2 

17. 

TE 

2-4 

11.7 

SS 

3-5 

15. 

Circular 

3-17 

11.6 

Circular 

4-7 

14. 

Axisymmetric 

4-12 

11.4 

Axisymmetric 

5-1 

14. 

Unsymmetric 

5-2 

17. 

Unsymmetric 

Assessment  and  evaluation.  From  the  above  fitting  and  prediction  exercises,  it  may  be 
concluded  that  the  inviscid  cap  model  generally  exhibits  good  fitting  and  predictive  capabili¬ 
ties  for  the  Colorado  concrete  data.  The  simulations  reported  herein  capture  the  overall  qual¬ 
itative  behavior  of  the  experimental  response. 
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COMPARISON  OF  EXPERI .  SIMUL.  DATA  FOR  CONCRETE  TEST 
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COMPARISON  OF  EXPERT.  &  SIMUL.  DATA  FOR  CONCRETE  TEST 
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SECTION  4 


CLOSURE 

A  systematic  estimation  procedure  for  the  parameters  involved  in  the  cap  model  to 
given  experimental  data  has  been  developed,  based  on  a  modified  Marquardt-Levenberg 
optimization  algorithm.  This  procedure  has  been  applied  to  the  extensive  experimental  pro¬ 
gram  carried  out  at  the  University  of  Colorado  and  reported  in  [4],  It  is  emphasized  that  due 
to  the  nonconveruional  character  of  this  experimental  data,  standard  fitting  procedures  (e.g., 
Desai  [5,6])  based  on  conventional  tests  can  not  be  employed.  The  numerical  simulations 
performed  on  the  basis  of  these  data  support  the  good  predictive  capabilities  of  the  cap  model 
for  concrete  materials. 
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APPENDIX  A 

LISTING  OF  PARAMETER  ESTIMATION  PROGRAM 


program  fit 

c  •*••*•*•••*••••« 

c....  Program  for  parameter  fitting  from  experimental  test  data 
c  for  the  cap  model 

implicit  double  precision(a-h.o-z) 
common/ fix/bulkm,shearm,zm 
common/aa/n(6),mm 
common/prop/Itype,tcut,fcut 
common  ystar<500,6),w(500,6) 
dimension  f(3000),q(21000),work(6063) 
c....  Specify  elastic  material  parameters  and  initial  cap  (z)  parameter: 
read(5,*)  bulkm,shearm,zm 
write(6,2004)  bulkm.shearm.zm 
c....  Input  size  of  experimental  steps 
read(5,*)  (n(j)j=l,6) 
mm  =  n(l)+n(2)  +  n(3)+n(4)+-n(5)  +  n(6) 
do  10j=l,6 
if(n(j)-gt.500)  stop  5 
10  continue 

c....  Input  Experimental  Data  for  y  (  sig-33) 
do  20  j=  1,6 

read(5,l000)  (ystar(lj),l=  l,n(j)) 

20  continue 

c....  Call  Optimization  Algorithm  ~ 
c  Modified  Levenberg-Marquardt  Algorithm, 
call  opt(f,q,work) 
stop 

c  Format  Statement 
1000  format(8fl0.0) 

2004  format(// 

*  tS, 'BULK  MODULUS  = ’,dl4.6/ 

*  t5, ’SHEAR  MODULUS  =  \dl4.6/ 

*  t5, ’INITIAL  CAP  (  Z  )  =  \dl4.6/) 
end 


c 


subroutine  opt(f.q.work) 


c 

c. 

c 


c.. 


..  Program  to  calc,  the  optimal  values  for  cap  model 
by  the  Modified  Levenberg-Marquardt  algorithm. 
..  The  optimization  criterion  is  with  respect  to  the 


c  least  square  norm. 

implicit  double  precision(a-h.o-z) 
external  func 

common/fix/bulkm,shearm,zm 

common/prop/ltype,tcut,fcut 

common/aa/n(6),mm 

common/bb/indi 

common/ha/old(7) 

common  ystar(500,6),w(500,6),y(500,6) 
dimension  parm(4),para(7),f(  1  ),q(mm,  1  ),g(28),work(  1 ) 
dimension  soss(6),sosorr(6),rmsso(6),rmsyy(6),reli(6) 
dimension  eigval(7),eigvec(7,7),wk(7) 
c..„  Parameters  for  IMSL  —  ZXSSQ 
ixjac-mm 

read(5,*)  nsig,eps,delta,maxfn,iopt 
iffiopt.eq.l)  read(5,*)  (parm(l), 1-1,4) 
c....  Initial  guess  for  parameters  alpha  "  w 
read(5,*)  (para(j)j=  1,7) 
write(6,2000)  (para(j)j=  1,7) 
c....  Input  weighting  matrix  W 
read(5,*)  iflag 
iftiflag.ne.  1)  then 
do  31  j- 1,6 
do  30  k=l,nG) 

30  w(kj)=l. 

31  continue 
else 

do  32  j=  1 ,6 

read(5,1000)  (w(kj),k=  l.nfj)) 

32  continue 
endif 

c....  Call  optimization  package  ZXSSQ 
indi  =  1 


c  (Save  the  old  parameters) 
do  40  k=  1,7 
40  old(k)  =  para(k) 


c 

call  zxssq(func.mm,7,nsig,eps.delta,maxfn,iopt, 

*  parm.para.ssq.f.q.ixjac.g.work.infer.ier) 
c....  Print  out  output  data 

writc(6,200l )  (para(k).k=  1.7) 
c.  ..  SOS  :  sum  of  squares  of  residuals 
c  SOSOR  :  sum  of  squares  of  observed  responses  (ystar) 
sos=0. 
sosor=0. 
do  50  i=  1.6 
sosorr(i)  =  0. 

50  soss(  i )  =  0. 
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do  100  k«l,n(j) 
diff»y(kj)-ystar(kj) 
soss(j)=soss(j)+diff**2 
sosorr(j ) = sosorr(j) + ystar(k  j  )**2 
100  continue 
sos-sos+soss(j) 
sosor  »  sosor + sosorr(j) 

200  continue 
c....  Overall  estimators 

rmssos  -  sqrt(sos/mm) 
rmsy = sqrt(sosor/ mm ) 
rel= rmssos/ rmsy 
c....  Test-j  estimators 
do  210  j- 1,6 
rmsso(j ) = sq  rt(soss(j  )/n(j )) 
rmsyy(j)=sqrt(sosorr(j)/n(j)) 
rell(j)=rmsso(j)/rmsyy(j) 

210  continue 
c 

write(6,2007)  sos 
write(6,2004)  rmssos 
write(6,2005)  rel 
c 

ml=n(l)+n(2) 
m2=ml+n(3) 
m3  =  m2  +  n(4) 
m4  =  m3  +  n(5) 
c....  Tesij 

do  220  j=  1,6 
write(6,2007)  soss(j) 
write(6,2004)  rmsso(j) 
write(6,2005)  rell(j) 

220  continue 
c 

c....  Compute  the  condition  number  for  G  =  =  Q  sup  T  *  Q  =  =  1/2  H 
e  Cond(G)  =  (max  lambda)  /  (min  lambda) 
call  cigrs(g.7,0.eigval.cigvcc.7,wk.icrr) 
if  (cigvaK  I  (.cq.O.OdO)  stop  'zero  cigval' 
cond=cigval(7)/cigval(  I ) 
wriie(6.3IUO)  cond 
return 
e 

1000  format(8fl  0.0) 

2000  format!// 

*  20.X.THE  INITIAL  GUESS  FOR  PARAMETERS’// 

*  IOX:aLPHA'.7X.THETA'.8X.'GAMA.8X.'BETA'. 

*  I  lX.'R'.l  IX.'D'.I  IX.'W'/ 

*  /5X.7FI2.6) 

2001  format)// 

*  20.\.'THE  OPTIMAL  VALUES  OF  PARAMETERS  ARE'// 

*  I0X.'ALPHA'.7X.'THETA'.8X. 'GAMA'. 8X. BETA'. 

*  I  lX.'R'.l  IX.'D'.I  IX.'W/ 

*  /5X.7FI  2.6) 


2004  format(//20x,’THE  TRUE  ROOT-MEAN-SQUARE  OF  PHI  =  \D15.7//) 

2005  FORMAT(//20x.’THE  NORMALIZED  RELATIVE  ERROR  -  \dl5.7) 
2007  format(///20x,’TRUE  SUM  OF  SQUARES  =  \dl5.7) 

3100  format(//20x.'CONDITION  NUMBER  OF  G  =  ’,dl5.7) 
end 


subroutine  func(para,in,ip,f) 
c  ******************************** 

c....  Function  evaluation  (stress  response)  and  residual 
c  computation. 

implicit  double  precision(a-h.o-z) 

common/fix/bulkm,shearm,zm 

common/prop/ltype,tcut,fcut 

common/aa/n(6),mm 

common/bb/indi 

c....  500:  max.  no.  of  data  pts  in  each  test 
c  6:6  strain  components 

c  6:6  tests 

common/ab/del(500,6,6) 

common/ha/old(7) 

common  ystar(500,6),w(500,6),y(500,6) 
dimension  para(  1  ),f(  1  ),delp(7),ytemp(  500),deltem(500,6) 
c....  Preserve  total  increments 
do  10  i- 1,7 

10  delp(i)-para(i)-old(i) 

c....  Check  if  constraints  are  violated: 
c  para(l)  >  0  required 
20  if  (para(  1  ).le.0.d0)  then 
go  to  30 

c  para(3)  >  0  required 

elseif  (para(3).le.0.d0)  then 
go  to  30 

c  para(3)  <  para(l)  required 

elseif  (para(3).gt.para(l))  then 
go  to  30 

c  para(3)  >  0. 1  *  para(l)  preferred 
elseif  (para(3).lt.0.  l*para(  1 ))  then 
go  to  30 

c  para(2)  >  0  required 

elseif  (para(2).!t.0.d0)  then 
go  to  30 

c  para(4)  >=  0.21  preferred 

elseif  (para(4).lt.0.21d0)  then 
go  to  30 

c  para(4)  <=  2  preferred 

elseif  (paru(4).gt.2.0d0)  then 
go  to  30 

c  para(5)>=  1.6  preferred 

elseif  (para(5).lt.  I  6d0)  then 
go  to  30 

e  para(6)  and  para(7)  >  0  required 

elseif  (para(6).le.0.d0.or.para(7).le.0.d0)  then 
go  to  30 


else 

c  ifO.K. 

go  to  50 
endif 

c....  Half  the  increments  for  parameters  if  constraints  are  violated. 
30  do  40  i=  1,7 
delp(i)=delp(i)/2. 

40  para(i)=old(i)+delp(i) 
go  to  20 

c....  Update  the  old  parameters 
50  do  60  i=  1,7 

60  old(i)=para(i) 
c 

do  70  j=  1,6 
if  (indi.eq.  1)  go  to  63 
do  62  k=  l,n(j) 
do  61  kk=  1,6 
deltem(k,kk)  *  del(k,kkj) 

61  continue 

62  continue 

63  call  main(ytemp,n(j),para,deltemj,indi) 
do  65  k=  l,n(j) 

y(kj)  =  ytemp(k) 
ytemp(k)=0.0 
do  64  kk=l,6 
del(k,kkj)  =*  deltem(k,kk) 
deltem(k,kk)  -  0.0 

64  continue 

65  continue 
70  continue 

indi  =  indi  + 1 
c 

knt-0 

do  200  j*  1,6 
do  100  i=  1  ,n(j) 
k=knt+i 

f(k)=(ystar(ij)-y(ij))*sqrt(w(ij)) 

100  continue 
knt»knt+n(j) 

200  continue 
return 
end 
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c  . . . 

subroutine  main(y,n,para,deI,ino,ind) 
c  ****•*•*•*******•****•*****•**•*••***•***********• 

c....  DEL  :  the  specified  strain  increment  vectors, 
implicit  double  precision(a-h.o-z) 
common/state/sig0(6) 
common/sta/geop,xint 
c....  Y  :  the  response  vector 
c....  PARA  :  the  parameter  vector 
dimension  del(500,6),sig(6) 
dimension  y(l),para(l) 
common/fix/bulkm,shearm,zm 
c....  Material  parameters  : 

common/prop/ltype,tcut,fcut 
common/elas/bulk, shear 
common/par  l/alpha,theta.gama, beta, r 
common/par2/d,w,z 

c....  Definition  for  parameters  ( just  for  convenience  ). 
bulk=bulkm 
shear=shearm 
alpha =para(l) 
theta =para(2) 
gama=para(3) 
beta=para(4) 
r=para(5) 
d=para(6) 
w=para(7) 
z  =  zm 

c....  IND  :  flag,  if  ,d=  1  ,  read  strain  increment  data 
c....  INO  :  identifier  for  test  #  ino  (  1-6  ). 
iffind.ne.  1)  go  to  200 
ifimo.ne.l)  go  to  50 

c....  Read  common  input  data:  ltype,tcut,sigO,geop,xint 
c...  Read  material  type  and  tension  cutoff  criterion 
c....  TCUT  is  in  terms  of  live  stresses. 
read(5,*)  ltype.tcut 

c....  Input  the  initial  states  of  stress  and  strain 
read(5,*)  (sig0(k),k=  1,6) 

c....  Input  the  geostatic  pressure  and  XINT(  the  initial  cap  ) 
read(5,*)  geop.xint 
c....  Strain  controlled  CAP  model 
c....  Read  input  data  del  and  initial  strain. 

50  dol00i=l,n 

read(5,*)  (del(i.k).k- 1,6) 

1 00  continue 

c....  Call  preprocessor  INITEL  to  calculate  elint  from 
c  given  xint  and  fcutfin  total  stress) 
c....  XINT  :=»  the  inital  X  value  for  the  initail  cap. 
c....  Z  :»  the  X  value  for  the  characteristic  initial  cap. 
c  i.e.  the  X  value  when  EVP  =  0. 

200  iffino.ne.  1)  go  to  210 

call  initel(xint, elint, nocon  1  ,nocon2) 
c....  Assign  initial  state  of  stress  and  strain  accordingly 


w 


subroutine  drv3D(n,y,del.el,sig) 


c....  This  routine  is  a  3-D  strain  history  driver  and  the 
c  variable  increments  are  deps  stored  in  array  del. 
implicit  double  precision(a-h.o-z) 
common/sta/geop 

dimension  del(500.6),sig(  1  ),deps(6),y(  1 ) 
common/prop/ltype,tcut,fcut 
common/elas/bulk,shear 
common/par  1  /alpha,theta,gama,beta,r 
common/par2/d,w,z 
c 

do  200  i-  l,n 
do  100  k- 1,6 
deps(k)«del(i,k) 

100  continue 

call  cap(sig,deps,geop,el,mtype,it,nocon,sj  1  ,sj2,xl,evpi 
*  ,ejl,ej2d,flejl) 
y(i)“sig(3) 


200  continue 
return 
end 
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subroutine  cap(sig,deps,geop.el,mtype,it.nocon,sj  I ,sj2, 
xl.evpi.ejl  ,ej2d,flej  1) 


c....  For  full  three-dimensional  stresses  and  strains 

c  computations  by  using  the  CAP  model 

c....  Strain  controlled  algorithm. 

c....  Stresses  and  strains  are  sig  and  eps,  respectively. 

c....  geop  =  geostatic  pressure  (overburden  stress) 

c....  el  =  hardening  parameter 

c....  mtype  :  0  *  tension  cutoff,  1  =  elastic,  2  =  failure, 

c....  3  =  cap  mode,  4=  cone  mode 

c....  it  »  #  of  iterations  for  CAP  mode  calculation 

c....  nocon  =  I  indicates  no  convergence  under  max  iterations 

c....  (nit)  restriction.  Otherwise  =  0 

c....  eps  =*  error  tolerance  parameter 

c....  ltype  :  1  =  soil,  2  =  rock 

implicit  double  precision(a-h,o-z) 

common/prop/ltype,tcut,fcut 

common/elas/bulk, shear 

common/par  l/alpha,theta,gama,beta,r 

common/par2/d,w,z 

dimension  sig(  1  ),deps(  1  ),s(6),de(6) 

data  eps/ 1  .d-6/ 

c.„.  Statement  function  for  exponential  with  negatively  large 
c....  argument  for  large  caps 

exps(z) =dexp(dmax  1  (-500., z)) 
c....  Failure  envelope  function  for  sj2 

fl(sj  l)=alpha-gama*exps(-beta*sj  l)+theta*sj  1 
d  1  (sj  1 ) = theta +gama*beta*exps(-beta*sj  l ) 
c....  Cap  statement  functions  for  f2  functional  forms 
c....  capl=»l(k) :  intersection  point  of  fl  &  f2, 
c....  x(k)  :  intersection  of  f2  &  j  1  axis 
capl(el)=dmaxl(0.0,el) 
ra(capi)=r 

x(el)=dmaxl(0.,el+ra(capl(el))*fl(el)) 

evp(xl)=w*(1.0-exps(d*(z-xl))) 

f2(sjl,xl,capi)  =  dsqrt((xl-capi)**2-(sjl-capi)**2) 

*  /ra(capi) 

c....  Elastic  moduli  functions 
bmod(sj  l,ev)=bulk 
smod(sj2,ev)*shear 


it=0 

nocon «0 

de  v = deps(  1 ) + deps(  2) + deps(  3 ) 

devb3=dev/3.0 

do  1  k=  1,3 

1  de(k)=deps(k)-devb3 
do  2  k  =  4,6 

2  de(k)=deps(k) 

press»(sig(  1  )+sig(2)+sig(3))/3.0 
do  3  k- 1,3 


do  4  k  =  4,6 

4  $(k)=sig(k) 
sjlt=3.*(press+geop) 
temp=0. 

do  11  k=  1,3 

1 1  temp=temp+0.5*s(k)*s(k) 
do  12  k=4,6 

12  temp  =  temp+s(k)*s(k) 
sj21=dsqn(temp) 

capi =capl(el) 
xl=x(cl) 
evpi  =  evp(xl) 

c....  Elastic  material  properties 
threek=3.*bmod(sj  lt.evpi) 
g=smod(sj21,evpi) 
twog=2.*g 
c....  Elatic  trial 

sj  1  =threek*dev+sj  1 1 
do  13  k=  1 ,6 

13  s(k)=s(k)+twog*de(k) 
ratio  =  1.0 

mtype=  1 

c....  Tension  limit  test 

tencut=dmaxl(fcut,tcut+3.*geop) 

if(sj  1  .gt.tencut)  go  to  10 

sjl=tencut 

ratio =0.0 

sj2=0. 

mtype=0 

c  If  no  contraction 

if(ltype.eq.2.or.el.le.0.0d0)  go  to  200 
c....  Tension  dilatancy  coding  for  soils  with  el.ge.0.0 
c....  Dilatancy  controlled  by  contracting  cap  up  to  el.ge.0.0 
ell=dmin  I  (0.0,el-eps*fl  (el)) 
xll=x(ell) 

denom =evp(xll)-evpi 
if(denom.lt.0.0d0)  go  to  5 
el=0.0 
go  to  200 

5  devp=dev-(sjl-sjlt)/threek 
denom =dmin  1  (denom, devp) 
el=el+devp*(ell-el)/denom 
el=dmaxl(0.0,el) 

go  to  200 

c....  Check  if  failure  envelope  mode  is  invoked 

10  continue 
temp=0. 
do  14  k=  1,3 

14  temp=temp+0.5*s(k)*s(k) 
do  1 5  k=4,6 

15  temp=temp+s(k)*s(k) 
c  Calc.  J2’E 

sj2=dsqa(temp) 

sj2e=sj2 
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c  If  cap  mode 

iftsj  1  .gt.capi)  go  to  40 
ej2d=sj2 

c....  TMISES  is  the  sj2  value  at  the  corner  point(tmises>-fj  1) 
tmises = f2(capi,xl,capi) 
ejl=sjl 
fjl -fl(sjl) 
flejl=fjl 

fe=sj2-dmin  1  (fj  1  .tmises) 
c  If  elastic 

iHfe.le.O.OdO)  go  to  200 

c....  If  k0<0  (small  cap)  ,  no  contraction  allowed, 
c  k=k0  and  J 1  =  J 1 E  (von  Mises  transition) 
if  (el.lt.0.0d0)  then 
mtype=*2 
go  to  30 
c....  For  k0>=0: 

c  If  JlE=L(kO),  Jl=L(kO)=JlE,  k=k0 
elseif  (dabs(sj l-capi).le.l.d-6)  then 
mtype=4 
go  to  30 
endif 

c....  Failure  envelope  surface  calculation  (  fl  ) 
mtype  =  2 
elold=el 

c....  Iterate  to  find  new  k  &  J 1  . 

call  proj(deps,el,sj l,2,sj2,nocon,it,threek,g) 
c....  Consistency  check  for  cap  model: 

if(ltype.eq.2.or.elold.eq.0.d0)  el=elold 
if(sj  1  .gt.el)  el  *  sj  1 

if  (Itype.eq.l.and.elold.gt.O.OdO)  then 
ifi(dabs(el-sjl).le.l.d-6)  mtype=4 
el  =  max(el.O.OdO) 
endif 

el  =  max(el.O.OdO) 

30  fjl=fl(sjl) 

sj2=dminl  (fjl,  tmises) 
ratio =sj2/sj2e 
go  to  200 

c....  CAP  mode  calculation 
40  if(sj  1  .gt.xl)  go  to  50 
c  If  elastic 

iffsj2.le.f2(sjl,xl,capi))  go  to  200 
50  mtype=3 

call  proj(deps, el, sj  1,3, sj2.nocon.it, threek.g) 

ratio=0.0 

if(sj2e.ne.0.0d0)  ratio = sj 2/sj  2e 
200  continue 
c....  Update  dev.  stresses  . 
do  300  k=  1.6 
300  s(k)=s(k)*ratio 
press=sjl/3.-geop 
c....  Calc,  live  stresses 
do  400  k=  1.3 
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400  sig(k)=s(k)+press 
do  410  k=4,6 
410  sig(k)  =  s(k) 

c....  calc.  X  and  vol.  plastic  strain  . 
xl=x(el) 
evpi  =  evp(xl) 
return 
end 

c  ******»*****«**********.***•******•***•**«*****•******•* 

subroutine  proj(deps,el,sjl.mtype,sj2,nocon,it,threek,g) 

(.  ******************************************************** 

c....  Subprogram  to  calc,  the  k  and  J1  iteratively  by  modified 
c  Regula  Falsi  Secant  Method. 

implicit  double  precision(a-h,o-z) 

common/prop/ltype,tcut,fcut 

common/elas/bulk, shear 

common/par  1/alpha, theta, gama, beta, r 

common/par2/d,w,z 

dimension  deps(l) 

data  nit/600/ 

data  eps/l.d-6/ 

c....  Statement  function  for  exponential  with  negatively  large 
c....  argument  for  large  caps 

exps(z) = dexp(dmax  1  (-500.  ,z)) 
c....  Failure  envelope  function  for  sj2 

fl  (sj  1  )=alpha-gama*exps(-beta*sj  1  )+theta*sj  1 
d  l(sj  1  )=theta+gama*beta*exps(-beta*sj  1) 
c.„.  Cap  statement  functions  for  f2  functional  forms 
c....  capl=l(k) :  intersection  point  of  f  l  &  f2, 
c....  x(k)  :  intersection  of  f2  &  jl  axis 
capl(el)=dmax  1  (0.0, el) 
ra(capi)=r 

x(el)  =  dmaxl(0.,el-i-ra(capl(el))*fl(el)) 
evp(xl)  =  w*(  1 .0-exps(d*(z-xl))) 

f2(sjl,xl,capi)=dsqrt((xl-capi)*(xl-capi)-(sji-capi)*(sj  1-capi)) 

*  /ra(capi) 

d2(sjl,xl,capi)=-(sj  l-capi)/ra(capi)/dsqrt((xl-capi)**2- 

*  (sjl-capi)**2) 

c....  Elastic  moduli  functions 
bmod(sjl,ev)=bulk 
smod(sj2,ev)=shear 

c*****«***«**. 

nocon=0 
it =0 
sjle=sjl 
sj2e=sj2 
xl  =  x(el) 
evpi=evp(xl) 

c....  Convergence  criterion 
conv=eps*0. 1 
c....  Failure  mode 
c  Initial  guess 

if  (mtype.eq.2)  then 
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xll = x(ell) 

dcvpl=evp(xll)-evpi 
sj  1 1 = sj  1  e-threek*devpl 
ql=-(sjll+l.)/(sjle+l.) 
c 

xlr=xl 

devpr=evp(xlr)-evpi 
sj  lr=sj  le-threek*devpr 
qr=(xlr-sjlr)/(xlr-jle) 
go  to  47 
c....  Cap  mode 
45  ell = el 
elr*sjle 

if(sj  le.ge.xl)  ql = (el-sj  1  e)/(el-xl) 

if(sj  1  e.lt.xl)  ql=2.*sj2e/(sj2e+f2(sj  le,xl,capi))-l  .0 

xr=x(elr) 

sj  lr=sj  le-threek*(evp(xr)-evpi) 
qr=(xr-sj  1  r)/(elr-xr) 

47  qold=0.0 

c....  Modified  Regula  Falsi  Method 
do  80  it*  1, nit 
c  Secant  method 

el=(qr*ell-ql*elr)/(qr-ql) 

xl=x(el) 

devp=evp(xl)-evpi 
sj  1  =sj  1  e-threek*devp 
capi=capl(el) 
iffmtype.eq.3)  go  to  48 
c  If  Failure  mode 
c  el  >=-l 

if(sj  1  .gt.el)  qc=-(sj  1  + 1  .)/(el-t- 1 .) 
if(sj l.le.sj le)  qc=(xl-sj  1  )/(xl-sj  1  e) 
if(sj  1  .gt.el.or.sj  1  .le.sj  1  e)  go  to  60 
sj2=fl(sjl) 
go  to  49 

c  If  Cap  mode 

48  continue 

if(sj l.ge.xl)  qc=(el-sj l)/(el-xl) 
if(sj  1  .le.capi)  qc=(xl-sj  1  )/(capi-xl) 
if(sj  1  .ge.xl.or.sj  1  .le.capi)  go  to  60 
sj2=*f2(sjl,xl,capi) 

49  if  (mtype.eq.2)  then 

c****  If  cone  mode  (at  corner  pt.)  inside  failure  mode1 
iffdabs(el-sjl).le.l.d-6)  then 
c....  Correct  treatment 

s!ope  =  (sj  l-sjle)/(sj2e-sj2)*g/(3.*threck) 
desp=devp/(3.*slope) 
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desp=devp/(3.*dl(sjl)) 

endif 

else 

desp=devp/(3.*d2(sj  1  ,xl,capi)) 
endif 

a=sj2-g*desp 
error=sj2e-a 
qc=error/(sj2e+a) 
c....  Convergence  criteria 

if(dabs(error).le.conv)  go  to  90 
60  if(qc.gt.0.0d0)  go  to  70 
if  (mtype.eq.3)  then 
c  k  too  large 
elr=el 
qr=qc 

iffqold.lt.O.OdO)  ql  =  0.5*ql 
else 

c  k  too  small 
ell «  el 
ql=qc 

iflqold.lt.O.OdO)  qr=0.5*qr 
endif 
go  to  80 


70  if  (mtype.eq.3)  then 
k  too  small 
ell=el 
ql=qc 

if(qold.gt.0.0d0)  qr=0.5*qr 
else 

k  too  large 
elr-el 
qr-qc 

iffqold.gt.O.OdO)  ql=0.5*ql 
endif 

80  qold=*qc 


c....  If  no  convergence  within  NIT  iterations: 
nocon* 1 

c  If  cap  mode: 

if  (mtype.eq.3)  then 
sjl=dminl(sjl,xl) 
if(sj  l .lt.capl(elr))  sjl  =capi 
sj2 =dmin  I(sj2e,f2(sj  1  ,xl,capi)) 
c  If  failure  envelope  mode: 

else 

sj  1  =dmin  1  (sj  1  ,el) 
endif 


c  89  continue 
90  return 
end 


subroutine  initel(xint,elint,nocon  1  ,nocon2) 
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c....  This  routine  uses  secant  method  to  find  initial 
c  value  of  el(hardening  parameter)  for  a  given 
c  initial  x(el)  value. 

c....  Also,  it  solves  FCUT.the  intersection  of  FI  and 
c  J 1  -axis. 

implicit  double  precision(a-h.o-z) 

common/prop/ltype,tcut,fcut 

common/elas/bulk,shear 

common/par  l/alpha,theta.gama,beta.r 

common/par2/d,w,z 

data  eps,nit/l.d-6,60/ 

c....  Statement  function  for  exponential  with  negatively 
c  large  argument  for  large  caps 
exps(z)=dexp(dmaxl(-500.,z)) 
c....  Failure  envelope  function  for  sj2 

f  1  (sj  1  )=alpha-gama*exps(-beta*sj  l)+-theta*sj  1 
c....  Cap  statement  functions 
capl(el)*dmax  1(0.0, el) 
ra(capi)=r 

x(el)*el+ra(capl(el))*fl(el) 
c....  Elastic  modulus  function 
bmod(sj  l,ev)=bulk 
c....  Find  initial  el 

c  Solve  f(k)  =  x(k)-xint=0  ,  not  related  to  Z. 
c....  xint  is  reset  so  that  within  the  convergence  criteria 
c  xint  is  positive  (  because  we  assume  x>0,  l(k)>=0.  ) 
c....  noconl  =  1  means  no  convergence  for  initial  el  iteration 
c  nocon2=  1  means  no  convergence  for  fcut  iteration 
c 

xint =dmax  1  (xint,eps*0.000 1  *bmod(0.,0.)) 
c....  Make  initial  guess  kO 
el0=xint*0.1 
Q0=x(el0)-xint 

c....  Make  second  initial  guess  kl 

el  1  =(xint-0. 1  *dmax  1  (dabs(xint).f  1  (xint)))*0.05 
c....  Set  up  convergence  criterion 
conv=dmin  1(  1  ,d-7,xint) 
c....  Secant  iteration 
do  100  it*  1, nit 
111  =x(ell)-xint 

if(dabs(fll).lt.conv.or.dabs(ell-elO).lt.conv)  go  to  200 

el2 = el  1  -fl  1  *(el  1  -el0)/(fl  1  -DO) 

el0=ell 

no=fii 

ell  =el2 
100  continue 
noconl  =  1 
200  elint =el  1 
c....  Find  fcut 
c  Solve  fl(fcut)=0 
c....  Make  first  initial  guess  for  fcut 
fcut=dmin  1(0., elint) 
dei*fl(fcut) 
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if(del.eq.0.d0)  go  to  600 
c....  Make  two  better  initial  guesses  for  fcut 
do  300  it»  l.nit 
elO=fcut-del 

n0=fl(el0) 

if(fl0.1t.0.d0)  go  to  400 
del=  10.*del 
fcut»el0 
300  continue 
c....  Secant  iterations 
400  do  500  it=  l,nit 
fl  1  =  f  1  (fcut) 

if(dabs(fll).lt.conv.or.dabs(fcut-elO).lt.conv)  go  to  600 
el2  -fcut-fl  I  *(fcut-el0)/(fl  1  -flO) 
el0=fcut 
fl0=(ll 
fcut=el2 
500  continue 
nocon2= 1 
600  return 
end 
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subroutine  dprint(y,nl,n2,name) 


Program  for  printing  response  y  (sig-33) 
implicit  double  precision(a-h,o-z) 
dimension  y(  1 ) 
character*6  name 

write(6,2000)  name 
c....  Print  out  8  columns  each  time. 

do  I00j*nl,n2,8 
c....  JH  :  the  right-most  index. 
jh-j+7 

if(jh.gt.n2)  jh=n2 
write(6,200 1 )  (n.n^jjh) 
write(6,2002)  (y(k),k=j  jh) 

100  continue 
return 
c....  Format 

2000  format(////20x,a6/ 

*  20x,’ 

2001  format(//8x,8il5) 

2002  format(/8x,8dl  5.7) 

end 
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INPUT  INITIAL  GUESS  FOR  MATERIAL  PARAMETERS: 

3.2  .09  1.0  .49  4.  0.004  0.2 

INPUT  OPTION  FOR  WEIGHTING  MATRIX  (0  :  WEIGHT  =  IDENTITY): 

0 

INPUT  OPTION  FOR  SOIL  (1)  OR  ROCK  (2);  AS  WELL  AS  TENSION  CUTOFF  VALUE 
I  -0.3 

INPUT  INITIAL  STRESS  STATE: 

0.  0.  0.  0.  0.  0. 

INPUT  GEOSTATIC  OVERBURDEN  PRESSURE  AND  INITIAL  CAP  POSITION: 

0.  16. 

INPUT  STRAIN  HISTORY  OF  6  TESTS: 

NO.  1 

-0.0000020  -0.0000587  -0.0000635  0.  0.  0. 

0.0002395  0.0002086  0.0001715  0.  0.  0. 

0.0004127  0.0004370  0.0004799  0.  0.  0. 

0.0003060  0.0003555  0.0003940  0.  0.  0. 

-0.0001604  -0.0002206  -0.0001891  0.  0.  0. 

-0.0003913  -0.0003993  -0.0003953  0.  0.  0. 

-0.0001187  -0.0001047  -0.0001095  0.  0.  0. 

-0.0000720  -0.0000330  -0.0000827  0.  0.  0. 

0.0001420  0.0000992  0.0000880  0.  0.  0. 

0.0004722  0.0004828  0.0004767  0.  0.  0. 

0.0005390  0.0006185  0.0006685  0.  0.  0. 

0.0004835  0.0005567  0.0005853  0.  0.  0. 

-0.0001502  -0.0002710  -0.0001843  0.  0.  0. 

-0.000 1722  -0.0002279  -0.0002293  0.  0.  0. 

-0.0001348  -0.0001654  -0.0001603  0.  0.  0. 

-0.0004710  -0.0004949  -0.0005204  0.  0.  0. 

-0.0001411  -0.0001104  -0.0001218  0.  0.  0. 

-0.0001776  -0.0001480  -0.0002248  0.  0.  0. 

0.0004733  0.0004259  0.0004466  0.  0.  0. 

0.0004973  0.0005316  0.0005792  0.  0.  0. 

0.0004338  0.0005852  0.0005394  0.  0.  0. 

0.0011258  0.0011877  0.0012631  0.  0.  0. 

0.0000354  -0.0000029  0.0003382  0.  0.  0. 

0.0000676  0.0000314  0.0002335  0.  0.  0. 

-0.0000763  -0.0000882  0.0001768  0.  0.  0. 

0.0000006  -0.0000269  0.0001764  0.  0.  0. 

0.0000336  0.0000646  -0.0000060  0.  0.  0. 

0.0000244  0.0000372  -0.0000439  0.  0.  0. 

0.0000506  0.0000774  -0.0000951  0.  0.  0. 
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0. 

0.0002324 

0. 

0. 

0. 

0.0001151 

0. 

0. 

0. 

0.0001662 

0. 

0. 

0. 

0.0003093 

0. 

0. 

0. 

0.0001864 

0. 

0. 

0. 

0.0001103 

0. 

0. 

0. 

0.0013489 

0. 

0. 

0. 

0.0017758 

0. 

0. 

0. 

BULK  MODULUS  =  0.2l0000d+04 
SHEAR  MODULUS  =  0.l70000d+04 
INITIAL  CAP  (  Z  )  =  O.OOOOOOd+OO 


THE  INITIAL  GUESS  FOR  PARAMETERS: 


ALPHA 

THETA 

GAMA 

BETA 

3.200000 

0.090000 

1.000000 

0.490000 

R 

D  W 

4.000000 

0.004000 

0.200000 

THE  OPTIMAL  VALUES  OF  PARAMETERS  ARE: 


ALPHA 

THETA 

GAMA 

BETA 

3.865751 

0. 100000 

1.163779 

0.443505 

R  D 

W 

4.433298 

0.003223 

0.429271 

TRUE  SUM  OF  SQUARES  =  0.6758328d+03 

THE  TRUE  ROOT-MEAN-SQUARE  OF  PHI  =  0. 1 537222d+0 1 

THE  NORMALIZED  RELATIVE  ERROR  =  0.2221052d+00 


CONDITION  NUMBER  OF  G  =  0.161 14l3d+05 
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